
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

'''
# [0,1) 均匀分布取样
x = np.random.random(10000)
#print(x)
#plt.hist(x, bins=np.arange(0,1,0.1), rwidth=0.8 )
#plt.show()
'''



'''
# 标准正态分布取样
help(np.random.randn)
x = np.random.randn(10000)
plt.hist(x, bins=np.arange(-5,5,0.1), rwidth=0.8 )
plt.show()
'''

# 任意分布取样
def f(x):
    return 1.0/np.sqrt(2*np.pi) * np.exp(-x*x/2)

#print( np.sqrt(quad(f, -100, 100)) ); exit(1)

# 返还 [a,b] 区间上服从分布 f(x) 的伪随机数
# fmax >= Max f([a,b])
def rand(f, a,b,n,fmax):
    rand1 = np.random.random(n)
    rand2 = a + (b-a) * np.random.random(n)
    return [ rand2[i] for i in range(n) if fmax * rand1[i]<f(rand2[i]) ]

plt.hist(rand(f,-3,3,100000,1), bins = np.arange(-3,3,0.1), rwidth=0.8, density=True)
# density=True 即绘制统计概率密度
x = np.arange(-3,3,0.1); y = f(x)
plt.plot(x,y)
plt.savefig("../../../doc/marp-slides/homemade.standard.distribution.png")
plt.show()

'''
# [-3,3]均匀分布中取30个数，再做平均,　理论上，这个平均值服从 N(0, 0.1)
x = [np.average(-3+6*np.random.random(30)) for x in range(10000)]
plt.hist(x, bins = np.arange(-1,1,0.03), rwidth=0.8, density=True)
x = np.arange(-1,1,0.03); y = 1.0/np.sqrt(2*np.pi*0.1) * np.exp(-x*x/2/0.1)
plt.plot(x,y)
plt.show()
'''